% rundtdS2.m
% 

%% clear and close if memory, display are issues
% clear all;
% close all;
%% options
load F18b
nfr = 3; nslice = 4;

%% setup
tic
recon.Ncpus = 1; % Number of CPUs available for parallel computing, '1' for serial computing

% characterize imaging experiment

pseq.fovcm = 22;  % field of view in cm

trajectory_file = 'KRT090320ab';NC=177;%circ 1.5x hexrosette 3.4g/cm,50ms,tweak 4.0,
%trajectory_file = 'KRT090401b';NC=177;

eval(['load ' trajectory_file]);  % trajectory data and eddy current phases
NL=floor(length(kss)/NC);
pseq.sw = 200000;  % sampling frequency
fid.delt=1/pseq.sw;
regfacm=0;regfacf=0;
ffac=100;
reso=88;
recon.alphaM = regfacm;  % regularization parameter for magnitude
recon.alphaf = regfacf;
%A
%  recon.NLIST = NC*[2;4;6;8;10;14;22;40;NL];
%  recon.FLIST = [.09;.07;.05;.03;.03;.02;.02;.01;.01];
%  recon.NLIST = NC*[2;4;6;8;10;14;22;40;NL];
%  recon.FLIST = [.07;.05;.04;.03;.03;.02;.02;.01;.01];
%B:
%  recon.NLIST = NC*[8;10;14;22;40;NL];
%  recon.FLIST = [.03;.03;.02;.02;.02;.01];
%  
recon.NLIST = NC*[8;10;14;22;40;NL];
recon.FLIST = [.04;.03;.03;.03;.03;.01];
wfi = 1;  % relative scaling of exponential imaginary
wfr = 1;  % relative scaling of exponential real
%=====================

model = model_setup(reso,pseq.fovcm,wfr+1i*wfi);
model.y = model.y + 2.0;%shift effective fov by 2 cm.
%model = model_setup_hex(reso,pseq.fovcm,wfr+1i*wfi);
recon.wfi = model.wfi;  % preconditioning weight for frequency (?)
recon.wfr = model.wfr;  % preconditioning weight for decay (?)
recon.NLENGTHS = size(recon.NLIST,1);  % # of outer loops in PLCG
%recon.FLIST = .02*ones(recon.NLENGTHS,1);
recon.NIT = 120;  % # of iterations
recon.GSS = 0;  % use golden section for line search
recon.del = 4e-09;  % standard step size at which to measure gradient in Cost function CC

recon.bounddecay = 5000/pseq.sw;  % threshold to bound decay parameter
recon.magdown = 0.2;  % reduce mag by this factor if decay exceeded
recon.plot = 0;  % show plots on the fly every recon.plot iterations (0=off)
recon.vis = 2;  % show images on the fly every recon.vis iterations (0=off)
recon.val = 0;  % show criterion values on the fly
recon.time = 0;  % show times on the fly
recon.mat = 'matrix';  % matrix, separable, or filter; save computation or memory

pseq.FR1 = real(squeeze(FID9218(241:11602,nslice,1,nfr)));
pseq.FI1 = imag(squeeze(FID9218(241:11602,nslice,1,nfr)));
pseq.FR2 = real(squeeze(FID9218(241:11602,nslice,2,nfr)));
pseq.FI2 = imag(squeeze(FID9218(241:11602,nslice,2,nfr)));
% phip and phir computed from calibration data using krdpenu19
pseq.krd = real(kss);  % readout k-coordinate
pseq.kpe = imag(kss);  % phase-encode k-coordinate
pseq.phir = phi;  % measured eddy current phase
pseq.phip = 0;  % measured eddy current phase 
pseq.dcoff = 0;
pseq.offr = 0;
kr=abs(kss);

%fid = fid_setupOO(pseq,recon.Ncpus,ffac,nfr);
fid = fid2_setupHMO(pseq,recon.Ncpus,ffac,nfr);
%fid = fid2_setupHNO(pseq,recon.Ncpus,ffac,nfr);
%fid = fid2_setupHMHM(pseq,recon.Ncpus,ffac,nfr);
model.pvec = zeros(1,3*length(model.ix),'single');


kbasis = make_kbasis_2(fid,model,recon.mat,pseq.fovcm,recon.Ncpus);
toc


%% run
whos
[model,kbasis] = cgparse6S_2(fid,recon,model,pseq,kbasis);

%% configure, save and plot results
modelsav = model;
model = vec2im_model2(model);
R2e = real(model.exp)*(-model.wfr/fid.delt);
frmap = imag(model.exp)*(model.wfi/(2*pi*fid.delt));
M0C1 = model.amp1;
M0C2 = model.amp2;
model = im2vec_model2(model);
